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Abstract 


This  paper  analyses  two  methods  based  on  the  Lanczos  algorithm  for  solving  large  sparse 
symmetric  linear  systems  with  several  right  hand  sides.  The  methods  examined  are  suitable  for 
the  case  where  the  right  sides  are  not  too  different  from  one  another  as  is  often  the  case  in  time 
dependent  or  parameter  dependent  problems.  We  will  show  in  particular  that  a  modified  Lanczos 
algorithm,  introduced  by  Pariett  is  in  some  sense  equivalent  to  the  block  Lanczos  algorithm. 
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1.  Introduction 

In  many  applications  one  needs  to  solve  several  symmetric  linear  systems  of  the  form 

AxW-b®  i=l,2..k.  (1) 

When  all  of  the  right  hand  sides  are  available  simultaneously,  then  several  block  methods  can  be 
successfully  applied  to  (1),  in  particular  the  block  Lanczos  algorithm  [5,  11],  the  block  Stiefel 
method  [11,  12]. 

In  practice  it  is  often  the  case  that  the  ’s  are  not  available  at  the  same  time,  i.e.  that  a 
given  right  hand  side  depends  on  the  solutions  x^,j=l,..i-l  of  the  previous  systems.  Then 
the  block  methods  are  no  longer  applicable.  For  this  common  situation,  Pariett  [7]  suggested  to 
use  the  Lanczos  algorithm  to  solve  the  first  system  and  to  save  the  Lanczos  vectors  thus 
generated  in  order  to  provide  good  approximate  solutions  to  the  subsequent  systems.  For 
example,  an  approximate  solution  to  the  second  linear  system  can  be  obtained  by  using  a 
projection  (Galerkin)  technique  onto  the  Kiylov  subspace  generated  when  solving  the  first  linear 
system.  We  refer  to  this  as  the  Lanczos-Galerkin  projection  procedure.  The  approximation 
obtained  from  the  Lanczos-Galerkin  projection  process  alone  niay  not  be  sufficiently  accurate  and 
a  further  refinement  is  often  needed.  A  suitable  and  efficient  way  of  improving  the  Lanczos- 
Galerkin  approximate  solution  is  to  use  a  special  Lanczos  process  introduced  by  Pariett  [7]  which 
consists  of  orthogonalizing  the  current  Lanczos  vector  not  only  against  the  previous  two  vectors, 
as  is  classically  done,  but  also  against  the  Lanczos  vector  of  the  previous  Krylov  subspace.  We 
will  refer  to  this  technique  as  the  modified  Lanczos  algorithm. 

The  purpose  of  the  present  paper  is  to  analyse  these  techniques,  from  the  theoretical  point  of 
view.  We  will  establish  an  error  bound  which  will  show  that  the  Lanczos-Galerkin  procedure  will 
provide  a  good  accuracy  under  the  condition  that  the  residual  vector  of  the  new  system  is  almost 
contained  in  the  previously  generated  Krylov  subspace.  We  will  also  show  that  the  modified 
Lanczos  algorithm  is,  in  a  certain  sense,  equivalent  to  a  block-Lanczos  method  with  a  particular 
initial  system. 

We  will  start  by  a  brief  presentation  of  the  Galerkin  projection  methods  based  on  the  Lanczos 
algorithm  as  described  in  [7].  Then  we  will  analyse  the  techniques  from  a  theoretical  point  of 
view  and  give  a  priori  error  bounds  for  the  Lanczos-Galerkin  method.  Finally  we  wilt  show  the 
relation  with  the  block  Lanczos  method. 
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2.  The  Lanczos  algorithm  for  solving  a  linear  system 

In  this  section  we  briefly  describe  the  Lanesos  method  for  solving  linear  systems.  Consider  the 
(single)  linear  system: 

Ax»b  (2) 

Suppose  that  a  guess  xQ  to  the  solution  is  available  and  let  rQ  be  its  residual  vector  rQ  =  b— AxQ. 
Then  the  Lanesos  algorithm  for  solving  (2)  can  be  described  as  follows: 


Algorithm: 

Stage  1:  Generate  the  Lanezos  Vectors 

•  Start:  Vj=r0/(/?:=||r0||) 

•  For  j=«=I,2..m  do: 


Compute: 


•i-(Ayi,Tj) 

(3) 

Vi A,i  -  "M’H  ' 

(  Vo  **0  ) 

(4) 

V.  HI  v.ll 

(5) 

Vi V  Ah 

(6) 

Stage  2  :  Form  the  approximate  solution 

*o  +  v„  T„'  <*,)  (7) 

where  Vm  ■■  [vj,v2fVm]  and  Tm  is  the  tridiagonal  matrix: 


°1  ^2 
°2  * 


(8) 


In  theory,  the  vectors  v.  computed  from  this  process  form  an  orthonormal  basis  of  the  Krylov 
subspace  Kffl  •—  span{r0,Ar0,..Am',r0}.  If  we  denote  by  Vffl  the  matrix  V  —{vlv.v  j  then  we 
have  V^AV  ««Tm  which  means  that  Tm  is  nothing  but  the  matrix  representation  of  the  section 
of  A  in  the  Krylov  subpace  Km  with  respect  to  the  basis  V  .  Furthermore,  it  is  easily  seen  that 
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the  Lane  Eos  algorithm  realizes  a  projection  process,  i.e.  a  Galerkin  process,  onto  the  Krylov 
subspace  Km,  see  e.g.  [7,  9].  In  other  words  the  approximate  solution  xm  can  be  found  by 
expressing  that  it  belongs  to  the  affine  subspace  x0+Km  and  that  its  residual  vector  b-Axm  is 
orthogonal  to  Km.  Denoting  by  Pm  the  orthogonal  projector  onto  Km,  this  means  that  the 
Lanczos  method  solves  the  approximate  problem: 

Find  x  €  xA+K  such  that: 
cd  u  m 

Pn,<b-  A*m>-0  (») 

The  approximation  thus  computed  is  identical  with  that  provided  by  m  steps  of  the  conjugate 
gradient  (CG)  method  when  A  is  positive  definite  [7].  When  A  is  not  positive  definite  this 
relation  between  the  Lanczos  algorithm  and  the  CG  method  can  be  exploited  to  derive  stable 
generalizations  of  the  CG  algorithm  to  symmetric  indefinite  systems  [6,  7,  2,  10]. 

A  well  known  and  troublesome  misbehavior  of  the  Lanczos  algorithm  is  the  loss  of 
orthogonality  of  the  Vj’s.  Fortunately,  this  does  not  prevent  the  method  from  converging  but 
often  results  in  an  important  slow  down.  Parlett  [7]  and  Simon  [13]  have  made  the  important 
observation  that  the  fast  convergence  properties  can  be  regained  by  resorting  to  different  sorts  of 
partial  reorthogonalizations.  This  important  matter  will  be  further  discussed  in  the  last  section. 

3.  The  Lanezos-Galerkin  projection  method 

Consider  now  the  two  linear  systems 

A  x(i)  —  b(i>  i— 1,2  (10) 

and  assume  that  m  steps  of  the  Lanczos  algorithm  described  in  the  previous  section  have  been 
performed  to  solve  the  first  system  in  a  first  pass.  We  would  like  to  use  the  information  gathered 
during  the  solution  of  the  first  system  to  provide  an  approximation  to  the  second  system: 

Ax,J,-b,J|.  (11) 

Clearly,  we  assume  that  the  vectors  Vj,  i— *1 ,2...m  as  well  as  the  tridiagonal  matrix  (8)  have  been 
saved,  possibly  in  some  secondary  storage. 

Suppose  that  we  have  an  initial  guess  x^  to  the  solution  of  (11).  Let  r^  be  the  residual  vector 
of  x(J,  that  is,  r^  -»  b^^-Ax^.  A  natural  way  of  improving  the  approximation  x^  is  by  means 
of  the  Galerkin  method  onto  the  Krylov  subspace  Km  generated  for  the  solution  of  the  first 
system.  This  will  yield  an  approximate  solution  s  defined  by: 
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«  -  *'?  +  vraTfnl  Vj  rlj  «  x(»j  +  y 
which  is  obtained  by  solving  the  m-dimensional  problem: 

Pm(  -  A  *  )  -  0 

or  equivalently 


(12) 


Vj(  b<2>  -  At)  =  0  (13) 

for  z  belonging  to  the  affine  subspace  x^  +Km-  We  will  refer  to  the  above  method  as  the 
Lanczos-Galerkin  process. 


One  important  question  concerning  the  approximation  obtained  from  the  above  Lanczos- 
Galerkin  process  is  its  accuracy.  Let  us  first  recall  that  when  A  is  positive  definite  we  can  define 
the  A”1- norm  of  a  vector  as  ||u||A-i— =(A_1x,x)1/2  and  that  the  projection  process  described  above 
minimizes  the  A-1- norm  of  the  residual  vector  over  all  vectors  of  the  form  x^  +  y  where  y 
belongs  to  Km  (4).  In  the  following  we  denote  by  Pm  the  orthogonal  projector  onto  Km,  where 
Km  is  the  Krylov  subspace  obtained  from  m  steps  of  the  Lanczos  algorithm  for  solving  the  first 
system. 

Proposition  1:  Assume  that  A  is  symmetric  positive  definite  with  largest  eigenvalue 
and  smallest  eigenvalue  X^.  Then  the  approximation  z  obtained  from  the  Lanczos 

Galerkin  projection  process  (12)  is  such  that 

l|b(!)  -A  i  llA-.  -  II  (I-P„)r|2j  llv.  +  <  04) 

where 


IM  - 

Tmb) 

in  which  7  —  A . 
[v ifAv.,..,Am'1Vj} 


first  ktnd. 


+  /  A/  *  a  is  the  first  component  of  Pmr^2J  in  the  basis 

and  Tm  represents  the  Chebyshev  polynomial  of  degree  m  of  the 


(15) 


Proof:  The  residual  vector  r  *■  b^  -  A  z  is  such  that 


F  —  b  -  A(x<2)  +  y) 

“  r*0  "  A  y 

where  y  is  the  vector  of  Km  computed  from  equation  (12).  The  residual  r  can  be  further 
decomposed  as: 

f-cv1?  -  Ay  )  +  (l-P„,)r<2>  (Id) 

Clearly,  the  projection  method  (12)  also  solves  the  system  A  j  ■  P„„r|(  by  the  same  Galerkin 


5 


process.  Hence  the  vector  y  also  minimizes  ||  Pfflr^  -Ay  ||A-i  over  all  vectors  y  of  Km.  Next, 
from  (16)  we  have  by  the  second  triangle  inequality: 

I II  r  ||A-,  -  llO-PJr1?  llA-«  I  <  II  -  AJ||A-,  (17) 

Let  m  Jet  <(?)  —  II  P„r<§'  -  Aj||a-i,  end  write  P_r*J»  at  —  qlA),,  where  q  is  a  polynomial 

of  degree  not  exceeding  m-1.  Clearly,  the  scalar  a  defined  in  the  proposition  is  the  constant  term 
in  the  polynomial  q,  i.e.  a=q(0).  Since  y  minimizes  r(y)  for  y  €  Kffl,  if  we  write  y=s(A)v1,  and 
define  the  polynomial  p(X)=q(X)-Xs(X),  we  see  that 


•  A  y||A-i 


min 

t  €  P*.,.  p(0)=» 


II  Vi^\\A-i 


(18) 


where  P  j  represents  the  set  of  all  polynomials  p  of  degree  <m-l.  The  above  equality  can  be 
rewritten  as: 

ii  p.*?  - A  yiU-‘  - 1 a  I  pem,i“  11  p^NHa-* 

The  last  term  of  the  right  hand  side  is  a  classical  factor  in  the  theory  of  the  conjugate  gradient 
method  and  a  well-known  upper  bound  for  it  is  available  (e.g.  [2])  and  yields 

H  pmr(?  * A  yiU->  ^  |a|  TTrT  (19) 

with  "t  =■  (Xj+Xn)/(Xj-Xn).  Tffe  result  finally  follows  from  inequality  (19)  and  inequality 

(17)P 

Let  us  now  interpret  the  result  of  the  proposition.  Notice  that  if  r^  belongs  to  the  previous 
Krylov  subspace  Km,  then  the  term  j|  (I-Pm)r^  ||A_i  in  the  right  hand  side  of  (14)  vanishes.  The 
proposition  then  indicates  that  in  this  case  the  method  will  provide  a  good  accuracy  when  a  is 
not  too  large.  In  fact  the  accuracy  will  be  of  the  same  order  as  that  obtained  from  m  steps  of  the 
regular  conjugate  gradient  method.  Note  that  if  a— >0  then  the  term  <  is  zero.  As  a  consequence 
an  extreme  case  where  the  new  system  can  be  exactly  solved  by  the  application  of  the  projection 
process  would  be  when  the  two  following  conditions  are  fulfilled: 

•  (M’J  •‘V  “  0.  i  '•  •'?  €K„ 

•  and  a— *0,  i.e.  r^  has  no  component  in  vt; 

The  opposite  extreme  case  is  when  the  projection  process  leaves  the  starting  approximate 
solution  unchanged.  This  happens  when  Pm  r^  —  0,  i.e.  when  is  orthogonal  to  Km-  In 


this  case  a  =  0  and  the  proposition  yields  ||r  ||A-i=|j  r^llA-i  which  is  clearly  true  since  i  =  x^. 

More  realistic  situations  arising  in  practice  will  lie  somewhere  between  these  two  extremes.  For 
such  general  cases  the  proposition  shows  that  the  error  eonsits  in  a  ‘small’  part  <  and  a  ‘large’ 
part  ||  (I~Pm)r^2Q  ||A-i  •  The  ‘small’  part  is  usually  as  small  as  would  be  obtained  from  m 
‘average’  steps  of  the  conjugate  algorithm.  The  ‘large’  part  depends  essentially  on  the  new  system 
and  can  be  quite  large  as  compared  with  «.  Perhaps  the  most  interesting  and  useful  such 
situations  arise  in  time  dependent,  or  more  generally  parameter  dependent  poblems,  in  which  the 
right  hand  sides  b^,  i=l,2,..  change  very  little.  Then  the  system  can  be  expected  to  be  solved 
relatively  accurately  because  the  ‘large’  term  ||(I-Pm)r^||A-i  becomes  small. 

When  ||(I'Pm)  r^  ||A-i  is  large,  then  it  is  likely  that  the  error  in  the  A-,-norm  sense  cannot  be 
decreased  below  |]  (I-PJrty  ||A-/  by  the  projection  process  (12)  alone.  This  means  that  some 
additional  refinement  must  be  applied. 


4.  Refining  the  Lancios-Galerkin  approximation 

Let  us  start  by  summarizing  the  essential  of  the  two  stages  of  the  process  described  in  the 
previous  section. 

1.  We  have  solved  the  first  linear  system  Ax^  — b^  by  the  Lanezos  method  and  this  has 
provided  us  with  a  Krylov  subspace  Km  of  dimension  m,  an  orthonormai  basis  Vm=[vj,v2,..vm] 
of  that  subspace  and  a  tridiagonal  matrix  Tm,  representing  the  section  of  A  in  Km,  with  respect 
to  this  basis. 

2.  We  are  faced  with  a  new  system  Ax^=b^,  for  which  an  initial  guess  x^  is  available.  This 
approximation  is  improved  by  use  of  the  Lanezos-Galerkin  projection  process  which  yields  the 
approximation: 

l  mm  x(2)  +  V  T-1  VT  r^ 

0  m  in  rn  0 

whose  residual  vector  is  r  —  b  -  A  z. 

We  have  at  this  stage  shown  only  how  to  improve  an  initial  solution  x^  to  the  second  system 
given  the  data  gathered  from  the  first  system  by  use  of  a  Galerkin  projection  process  on  Km. 
The  accuracy  of  z  thus  obtained  may  be  far  from  sufficient  as  is  shown  by  the  comments  at  the 
end  of  the  previous  section.  We  are  therefore  faced  with  the  problem  of  improving  the  new 
approximation  z.  An  easy  way  of  achieving  this  is  by  a  completely  new  sequence  of  Lanezos 
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approximations  starting  with  z.  However,  this  does  not  fully  take  advantage  of  the  available 
information,  i.e.  of  the  previous  Krylov  subspace  and  the  representation  Tm.  A  suitable 
alternative  would  be  one  which  continues  naturally  the  previous  process  by  constructing  a 
sequence  of  subspaces  containing  K(n  and  of  increasing  dimension. 

One  such  process  was  first  introduced  by  Parlett  [7]  and  was  later  rediscovered  by  Carnoy  and 
Geradin  (l]  in  a  different  context.  The  following  algorithm,  which  will  be  referred  to  as  the 
modified  Lanczos  algorithm,  differs  only  in  its  presentation  from  Parlett’s  algorithm  and  the 
algorithm  in  \cite{Carnoy-Geradin}.  Its  purpose  is  to  compute  a  sequence  of  vectors  Wj,  i=l,2,... 
which  is  orthonormal  and  also  orthogonal  to  the  v.’s,i«"»l,..,m,  generated  for  the  first  system. 

Modified  Lanczos  Algorithm 


1. Start:  take  w(  —  r/||r1| 


2 .Iterate:  For  j— 1,2, ...do 
Compute  Oj=(Aw^Wj) 

*j=(Awj’Vn,) 

w]+1  -  Aw.  -  ajwj  -  ?.w KS.rm 

^j+l  ™  I^j+lH 

wj+i  “  wi+i/^j+i 


(20) 


All  the  difference  with  the  usual  Lancios  algorithm  is  that  at  each  step  we  now  orthogonalize 
against  one  more  vector,  namely  the  vector  v  .  We  claim  that  this  simple  modification  of  the 
Lancios  algorithm,  yields  a  sequence  of  vectors  so  that  the  system  Wp={v1,v2,....vm,w1,..wp}  is 
orthonormal.  This  is  an  important  property  since  it  will  allow  us  to  realize  in  a  simple  way,  as 
will  be  seen  shortly,  the  Galerkin  process  onto  the  span  of  Wp. 

Proposition  2:  Suppose  that  p  steps  of  the  Modified  Lanczos  algorithm  are  feasible. 

Then  the  sequence  forms  an  orthonomal  sequence  of  vectors. 

Proof:  Since  (vj}.— j  m  is  an  orthonormal  system,  we  must  prove  that  for  j=l,2,..  we  have: 

1.  Wj  is  orthogonal  to  the  v/s  i**l,.m 

2.  w.  is  orthogonal  to  the  previous  w/s,  i— l,2,..j-l. 
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The  proof  in  by  induction.  Clearly  the  above  property  is  true  for  j=l,  because  Wj  is  equal  to  r 
apart  from  a  multiplicative  constant  and  r  is  known  to  be  orthogonal  to  the  subspace  Km  by  the 
Galerkin  condition.  Suppose  that  the  propoerty  is  true  for  j  and  let  us  prove  that  it  is  true  for 
j+1,  i.e.  that 

(wj+l’Ti>  "  °-  i-l^~m  (21) 

(wj+1,wj)  -  0,  i=*l,2,..j.  (22) 

Consider  (21)  first.  By  construction,  (wj+i»Tm)“0»  80  we  can  restrict  ourselves  to  the  case  i<m. 

h+l^H^j+ir1  [  <Awi'Ti)  •  1 

By  the  induction  hypothesis  and  because  {▼;};_!  m  is  orthogonal,  the  terms  (w.,v.),  (w.  ^Vj)  and 
(v^,Vj)  on  the  right  hand  side  all  vanish.  The  remaining  term  (Aw.,v.)  can  be  expanded  as 
follows: 

(Aw.,vj)=(wj,Avi)  =  ( Vi+iTi+l+ttiTi+4vi-l) 

=-^+l(Wj,Ti+i)+ai(w.,Vi)+)Ji(Wj,Vi.i) 

Using  again  the  induction  hypothesis,  we  see. that  all  these  terms  are  in  turn  equal  to  zero.  This 
completes  the  proof  of  (21). 

Now  consider  (22). 

(Wj+l’WiH*  I  *j  +  l  l'1  I  (Awj>wi)  -  «j(wj(Wi)  -  Jf wH,Wj)  -«j(vm,wi) ) 

Assume  first  that  i<j  —  1.  By  the  induction  hypothesis  we  obviously  have 

(w.,w.)  =(w.  ,,w  )  =  (v  .w.)=  0. 

'  j’  v  '  j-1’  r  '  m’  i' 

Proceeding  as  before,  we  expand  the  remaining  term  (Aw.,w;)  as  follows: 

(Aw.,wi)=(wj,Aw.)  -  ( Vi+l^+l+aiwi+^iwi-I  +  fjvJ 
“  ^i+ 1( wj’wi+ 1  )+®i(wj  -wi)+^( i )  + 

which,  by  a  final  application  of  the  induction  hypothesis  shows  that  (AWj,Wj)=0.  Hence 
(Wj+i,Wj)=0,  for  i<j  -1.  For  i=j  and  i=j— 1,  the  scalars  a.,  and  ^  have  been  precisely  chosen 
so  that  the  property  is  true.  This  completes  the  proof.Q 


Consider  the  subspace  spanned  by  the  orthonormal  system  Wp=s[vj,v2,...vm,Wj,..wp)  which  we 
will  denote  by  Km  .  Note  that  Kffl  p  is  no  longer  a  Krylov  subspace  but,  as  will  be  seen  shortly, 
the  case  p»m  is  of  particular  importance.  The  matrix  representation  W^AWp  of  the  section  of 
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A  in  the  subspace  Km  ,  with  respect  to  the  basis  Wp  is  the  (m+p)  x  (m+p)  matrix: 


Q1  ^2 

a2  &3 
&3  a3 


■  0. 

0M 


"j  0_2  - 

02  -®2  03 

03  °3 


6 

—  -  I 


f  /* 

p  a 
p  P 


7' 


Hence  the  new  approximate  solution  obtained  at  step  p  of  the  projection  process  onto  the 
subspace  Km  is  given  by: 

,  -1  + WT-*  WjF 

p  p  m,p  p 

Noticing  that  W%=  wj  ||z||  Wj,  this  simplifies  into: 

*p  -  *  +  IMI  WP  T„-,'p 

where  «m+1  is  the  vector  of  length  m+p  whose  components  are  zero  except  the  (m+l)st  which  is 

equal  to  1.  Once  the  augmented  tridiagonal  system  Tm  s  =  em+i  ‘9  s°lve(i,  the  linear 

combination  W  s  of  the  vectors  v-  and  w.  should  be  formed  and  added  to  z.  The  whole  set  of 
P  *  ' 

vectors  Wp  must  therefore  be  kept  in  storage. 

It  is  important  to  interpret  the  method  outlined  above,  in  order  to  compare  its  rate  of 
convergence  with  that  of  the  regular  Lanczos  process.  We  will  establish  the  following  result  for 
the  particular  case  where  p=m. 

Proposition  3:  When  p=m,  the  modified  Lanezoa  process  ia  equivalent  to  m  steps  of 
the  block  Lanezoa  method  with  block  dimension  of  2,  with  starting  block  consisting  of 
the  vectors  v ;  and  wr 

Proofs  The  proof  amounts  to  showing  that  the  two  methods  realize  the  Galerkin  process  on 
the  same  subspace. 

For  the  block  Lanczos  method  [5,  14,  12],  the  subspace  is  simply  span{Sj,ASI,...Am‘lSj}  where 

s, 


The  modified  Lanczos  algorithm  is  a  projection  process  onto  the  subspace 


>an{vj,v2,...vm,Wj,w2....wm}.  ^  '*  a  si®ple  exercise  on  proofs  by  induction  to  show  from  the 
gorithm  that  the  Vj’s  and  w/s  are  linear  combinations  of  the  vectors  A^Vj  and  A^Wj, 
<k,j<m-l  and  vice-versa.  Thus  the  two  subspaces  are  identical  and  the  proof  is  complete.^ 

The  proposition  asserts  that  there  are  two  ways  of  realizing  the  block  Lanczos  method.  One  is 
le  regular  algorithm  in  either  its  block  form  [14]  or  its  banded  form  [8],  and  is  suitable  for  the 
ise  when  the  linear  systems  are  avaibale  simultneously.  The  other  way  is  the  modified  Lanczos 
Igorithm  which  is  more  suitable  when  the  right  hand  sides  are  not  available  simultaneously,  i.e. 
hen  the  right  hand  side  bl2)  depends  on  the  solution  x^  of  the  first  system. 

The  rate  of  convergence  of  the  block  Lanczos  algorithm  for  solving  linear  systems  has  been 
:udied  in  [5,  11]  and  we  will  not  report  the  results  here.  It  suffices  to  say  that,  not  surprisingly, 
lie  bounds  on  the  rate  of  convergence  of  the  block  method  are  superior  to  those  of  the  single 
ector  method.  We  should  point  out  however  that  our  experience  reveals  that  when  only  one 
ystem  is  to  be  solved  it  is  not  in  general  effective  to  use  a  block  method  as  a  means  of 
ccelerating  the  convergence  [11]. 

To  summarize,  the  modified  Lanczos  method  has  the  advantage  of  the  rapid  convergence  of 
be  block  Lanczos  method  without  the  drawback  of  requiring  that  the  second  right  hand  side  be 
vailable  at  the  same  time  as  the  first. 

.  Practical  Considerations 

One  important  feature  of  both  the  Lanczos-Galerkin  Process  and  the  modified  Lannczos 
rocess  is  that  We  must  save  a  large  number  of  vectors  in  some  secondary  storage.  This  may 
rem  impractical  at  first  but  there  are  numerous  reasons  why  it  is  not  always  so: 

•  Once  a  vector  has  been  computed  it  is  not  needed  until  the  convergence  of  the  process 
is  reached.  There  exists  a  simple  formula  for  evaluating  the  residual  norm  of  the 
solution  without  even  having  to  compute  the  solution  [7,  10]  thus  allowing  to 
determine  when  to  stop. 

•  There  are  supercomputer  systems  with  very  fast  auxiliary  memories,  e.g.  the  Cray- 
XMP  with  Solid-state  Storage  Device  (SSD). 

•  In  many  cases  the  dominant  cost  is  the  matrix  by  vector  product  and  therefore  the 
priority  is  to  economize  on  the  number  of  matrix  by  vector  multiplications.  The 
Lanczos-Galerkin  process  of  section  3  requires  only  one  matrix-vector  product  (for 
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computing  the  initial  residual  r^). 

The  Lanczos-Galerkin  process  was  successfully  used  in  the  context  of  stiff  ODC's  [3].  In  some 
ODE  problems  the  cost  of  a  matrix  by  vector  multiplication  can  be  extremely  high  and  the 
Lanczos-Galerkin  process  becomes  very  attractive. 

In  the  context  of  symmetric  generalized  eigenvalue  problems,  a  technique  similar  to  the  one 
described  in  the  previous  section  has  recently  been  presented  by  Carnoy  and  Geradin  [l],  who 
report  some  interesting  numerical  results. 

It  seems  important  that  the  vectors  v-’s  that  must  be  saved  from  the  first  linear  system,  be 
orthogonal  because  the  Lanczos-Galerkin  process  is  essentially  based  on  the  orthogonality  of  these 
vectors.  The  Selective  Orthogonalization  [7]  and  the  Partial  Orthogonalization  [13]  methods  can 
both  be  extended  to  the  modified  Lanezos  method  and  are  attractive  alternatives  to  the  more 
expensive  full  reorthogonalization  schemes.  Simon  [13]  has  shown  that  any  partial 
reorthogonalization  that  guaranties  semi-orthogonality,  i.e.  orthogonality  within  the  square  root 
of  the  machine  accuracy  t,  will  also  deliver  an  approximate  solution  vector  that  is  within  /T  of 
the  ideal  solution  vector  from  the  Krylov  subspace. 

Although  not  explicitely  mentioned,  it  is  clear  that  the  techniques  of  the  previous  two  sections 
can  easily  be  extended  to  systems  with  more  than  two  right  hand  sides.  The  resulting  algorithms 
are  straightforward  and  so  is  the  theory. 
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